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Abstract 

In  this  thesis,  we  develop  a  numerical  method  in  order  to  approximate  the  solutions 
of  one-dimensional,  non-linear  absorption-diffusion  equations.  We  test  our  method  for 
accuracy  against  a  linear  diffusion  equation  with  a  solution  that  can  be  written  in  closed 
form.  We  then  test  various  types  of  diffusion  and  absorption  terms  to  determine  which 
ones  produce  extinction  in  finite  time. 

We  also  develop  a  numerical  method  to  computationally  solve  diffusion-free  equa¬ 
tions.  We  compare  the  numerical  solutions  of  the  one-dimensional,  non-linear  absorption- 
diffusion  equation  and  the  diffusion-free  equation  and  we  find  that  for  the  cases  tested,  the 
numerical  absorption-diffusion  solutions  are  always  less  than  the  numerical  diffusion-free 
solutions.  Furthermore,  we  find  this  is  true  for  the  cases  tested  when  there  is  finite  and 
infinite  extinction  time. 

We  also  look  at  the  open  problem  where  we  have  slow  diffusion  and  weak  absorption 
but,  their  combined  effect  is  strong.  Our  results  provide  some  insight  into  the  answer  of 
this  problem. 
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FINITE  EXTINCTION  TIME  FOR  NON-LINEAR  ABSORPTION-DIFFUSION 


EQUATIONS 

/.  Introduction 

1.1  Motivation 

“Diffusion  is  the  process  by  which  matter  is  transported  from  one  part  of  a  system  to 
another  as  a  result  of  random  molecular  motions  [4] Despite  the  fact  that  molecules  have 
no  preferred  direction  of  motion,  there  is  an  overall  movement  of  molecules  from  regions 
of  higher  concentration  to  regions  of  lower  concentration.  This  can  be  explained  by  the 
fact  that  over  a  given  interval  of  time,  on  the  average,  a  definite  fraction  of  molecules  from 
the  higher  concentration  will  move  into  the  lower  concentration,  and  the  same  fraction 
of  molecules  will  move  from  the  lower  to  higher  concentration.  So,  simply  because  there 
are  more  molecules  in  the  higher  concentration,  there  is  a  net  transfer  of  molecules  from 
the  higher  to  the  lower  concentration,  as  a  result  of  random  molecular  motions  [4].  In 
addition  to  concentrations,  the  matter  moving  can  also  be  considered  as  heat  flow.  Thus, 
the  diffusion  equation  is  commonly  referred  to  as  the  heat  equation. 

The  n-dimensional,  linear  diffusion  equation  is  given  by 

ut{x,t)  =  a‘^'V‘^u{x,t),  X  G  0<t<oo,  (IT) 

where  U  is  usually  a  bounded  domain  in  A  linear  equation  means  that  the  dependent 
variable  u  and  all  its  derivatives  appear  in  a  linear  fashion,  i.e.,  they  are  not  multiplied 
together  or  squared  or  etc.  [6].  Equation  (1.1)  relates  the  rate  of  change  in  the  temperature 
or  concentration  profile,  u{x,  t),  with  respect  to  time,  uj  =  ^,  and  -|-  •  •  •  -|- 

1^,  which  in  one-dimension  is  the  concavity  of  the  temperature  or  concentration  profile. 

in  essence  compares  the  temperature  or  concentration  at  one  point  to  the  temperature 
or  concentration  at  neighboring  points  [6].  The  quantity  >  0  is  the  diffusion  coefficient 
measured  in  length^/time. 
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When  dealing  with  an  internal  heat  source,  an  internal  heat  sink,  or  an  irreversible 
reaction,  f{x,t,u{x,t)),  the  material  is  either  being  supplied  with  heat,  heat  is  being 
removed,  or  a  reaction  is  taking  place  at  the  location  x  and  instant  of  time  t,  and  may 
depend  on  the  temperature  or  concentration  profile,  u{x,  t).  We  call  /  the  absorption  term. 
This  leads  to  a  modification  of  the  n-dimensional,  linear  diffusion  equation  as  follows. 

ut{x,t)  =  a‘^'V‘^u{x,t)  —  f{x,t,u{x,t)),  x  G  Q,  0<t<oo,  (1-2) 

where  all  the  variables  are  the  same  as  in  equation  (1.1).  Equation  (1.2)  is  often  called  the 
absorption-diffusion  equation. 

Throughout  the  remainder  of  this  thesis,  we  will  consider  a  non-linear  absorption- 
diffusion  equation,  which  has  the  form 

ut  =  V^d>(u)  —  F{u),  X  G  fl,  0  <  t  <  oo.  (1-3) 

That  is,  the  absorption  term  depends  only  on  the  temperature  or  concentration.  We  will 
derive  the  one-dimensional  case  of  equation  (1.3)  in  the  next  section. 

1.2  Derivation 

We  will  consider  a  one-dimensional,  non-linear  absorption-diffusion  equation  with 
the  form 


ut{x,t)  =  ^{u{x,t))xx  —  F{u{x,t)),  0  <  X  <  L,  0  <  t  <  oo.  (1-4) 

We  will  formally  derive  equation  (1.4)  from  a  heat  flow  standpoint.  From  this  vantage, 
conservation  of  energy  is  our  cornerstone. 

Let  us  suppose  we  have  a  thin  rod  of  length  L  lying  centered  along  the  x-axis  from 
0  to  L.  We  assume  the  rod  is  thin  so  that  the  temperature  at  all  points  of  a  cross  section, 
A,  are  constant.  This  essentially  means  that  we  have  a  one-dimensional  rod.  Also,  let  us 
assume  that  the  rod  is  laterally  insulated  so  that  heat  can  only  flow  in  the  x-direction. 
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0 


X 


X  +  Ax 


L 


\ - 5^  X 

Figure  1.1.  Thin  Conducting  Rod. 

Dividing  the  rod  into  small  segments  of  length  Ax,  we  can  apply  the  principle  of 
conservation  of  heat  to  a  segment  of  the  rod,  [x,  x  +  Ax].  Therefore, 


Net  change  of  heat  inside  the  segment  [x,  x  +  Ax] 

=  Net  change  of  heat  across  the  boundaries  at  x  and  x  +  Ax  (1-5) 

+  Total  heat  generated  or  absorbed  inside  the  segment  [x,x  +  Ax]. 


Now,  at  any  time  t. 


Total  heat  inside  the  segment  [x,  x  +  Ax] 


nx-\-l\x 


Ac{u)p{u)u{s,  t)ds, 


(1.6) 


where  c{u)  is  the  thermal  capacity  of  the  rod  and  p{u)  is  the  density  of  the  rod,  both  of 
which  we  assume  depend  on  the  temperature.  Thermal  capacity  is  a  material  property  of 
the  rod  which  measures  its  ability  to  hold  heat. 

This  makes  the  change  of  heat  inside  the  segment  over  time  equal  to  the  derivative 
with  respect  to  t  of  equation  (1.6)  as  shown  below. 


m 


Net  change  of  heat  inside  the  segment  [x,  x  +  Ax] 

X-\-/\x  rX-\-/\x  Q 

Ac{u)p{u)u{s,t)ds  =  A  /  —[c{u)p{u)u{s,t)]ds. 

J  X 


(1.7) 


Next,  we  need  to  find  the  change  in  heat  across  the  boundaries  at  x  and  x  +  Ax. 


Net  change  of  heat  across  the  boundaries  at  x  and  x  +  Ax 
=  A{k{u{x  +  Ax,  t))ux{x  +  Ax,  t)  —  k{u{x,  t))ux{x,  t)}, 


(1.8) 
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where  k{u)  is  the  thermal  conductivity  of  the  rod  and  Ux{x,t)  is  the  thermal  gradient  at 
the  boundaries  at  x  and  x  +  Ax.  Thermal  conductivity  is  the  ability  of  the  material  to 
allow  heat  to  diffuse,  and  in  our  case,  is  dependent  on  the  temperature  of  the  material. 

Finally,  we  must  deal  with  the  heat  generated  or  absorbed  inside  the  segment. 


Total  heat  generated  or  absorbed  inside  the  segment  [x,  x  +  Ax] 

rx+Ax 


=  A 


f(u(s,t))ds, 


(1.9) 


where  f{u{s,t))  is  a  function  of  how  the  material  will  generate  or  absorb  heat. 
Substituting  equations  (1.7),  (1.8),  and  (1.9)  into  equation  (1.5)  we  get 

/•X+Ax  Q 

A  /  —[c{u)p{u)u{s,t)]ds  =  A[k{u{x  +  Ax,t))ux{x  +  Ax,t) 

J  X 

rx-\-Ax 

-  k{u{x,t))ux{x,t)]  +  A  /  f{u{s,t))ds. 

J  X 

(1.10) 

At  this  point,  we  will  assume  that  f{u{s,  t))  is  a  heat  sink,  i.e.,  a  source  of  heat  absorption 
in  the  segment,  instead  of  a  heat  generator.  Thus,  f{u{s,t))  <  0.  However,  we  will  change 
the  sign  in  front  of  f{u{s,t))  so  that  /  itself  is  nonnegative. 

Dividing  each  side  of  equation  (1.10)  by  A,  and  applying  the  mean  value  theorem  for 
integrals,  we  obtain 


d 

—  [c{u{i,t))p{u{£,,t))u{i,t)]Ax  =  k{u{xA  Ax,t))ux{x  + Ax,t) 

-  k{u{x,t))uxix,t)  -  f{u{^,t))Ax, 


(1.11) 


for  some  x  <  ^  <  x  +  Ax. 

Now,  dividing  both  sides  of  equation  (1.11)  by  Ax  gives  us 


^[c{u{^,t))p{u{C,t))u{^,t)]  = 


k{u{x  +  Ax,  t))ux{x  +  Ax,  t)  —  k{u{x,  t))ux{x,  t) 


Ax 


(1.12) 
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Letting  Ax  ^  0  in  equation  (1.12)  results  in 
d 

—  [c{u{x,t))p{u{x,t))u{x,t)]  =  [k{u{x,t))ux{x,t)]x  -  f{u{x,t)) 

or 

d 

—  [c{u)p{u)u\  =  [k{u)u^]a;  -  f{u). 

Next,  we  let  g{u)  =  k{z)dz  so  that  g{u)xx  =  [k{u)ux]x-  This  gives  us 

d 

—  [c{u)p{u)u]  =  g{u)xx  -  f{u). 


(1.13) 


(1.14) 


(1.15) 


Now,  we  let  w  =  c{u)p{u)u  and  we  assume  that  this  relationship  is  invertible  so  that  we 
can  solve  for  u  to  get  u  =  h{w)  for  some  function  h.  Substituting  this  into  equation  (1.15) 
we  obtain 

=  g{h{w))xx  -  f{h{w)).  (1.16) 

Letting  ^(w)  =  g{h{w))  and  F{w)  =  f{h{w))  in  (1.16)  we  arrive  at 


wt  =  ^>(w)a;a;  -  F{w), 


(1.17) 


our  one-dimensional,  non-linear  absorption-diffusion  equation. 

Finally,  we  switch  back  to  using  u  instead  of  w,  and  note  that  the  rod  lies  in  the 
interval  [0,  L]  to  get  0  <  x  <  L.  This  gives 

ut  =  ^{u)xx  —  T(rt),  0<x<L,  0<t<oo.  (1-18) 

The  analogous  problem  in  n-dimensions  is  given  by 

Ut  =  V‘^^{u)  —  F{u)  inflx(0,oo),  (1-19) 

where  is  a  bounded  domain  in  M"". 
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1.3  Problem  Statement 


Our  goal  for  this  thesis  is  to  determine  the  types  of  diffusion  terms,  ‘h,  and  absorption 
terms,  F,  for  equation  (1.19),  which  produce  extinction  in  finite  time.  We  would  also  like 
to  determine  which  conditions  on  and  F  do  not  produce  extinction  in  finite  time. 

Definition  1.3.1  Finite  extinetion  time  is  a  finite  time  to  >  0  sueh  that  u{x,t)  =  0  for 
all  {x,t)  G  n  X  [to,oo). 

Throughout  the  remainder  of  this  thesis,  we  will  assume  that  and  F  are  nonde¬ 
creasing,  nonnegative  C'^((0,oo))  n  C'([0,oo))  functions  that  satisfy 

$(0)  =  0,  F(0)  =  0,  ^>(s)  >  0,  F{s)  >0  if  s  >  0.  (1.20) 

A  (^([O,©©))  function  is  a  real-valued  function  that  is  continuous  on  the  interval  [0,oo). 
A  (^^((O,  oo))  function  is  a  real- valued  function  whose  first  derivative  is  continuous  on  the 
interval  (0,oo). 

The  assumptions  in  (1.20)  have  physical  significance.  ^(0)  =  0  means  that  when 
there  is  no  heat,  i.e.,  temperature  =  0,  there  is  no  diffusion.  T(0)  =  0  means  that  there  is 
no  absorption  if  there  is  no  heat.  d>(s)  >  0  and  F(s)  >0  for  s  >  0  means  that  diffusion 
and  absorption  occur,  respectively,  when  there  is  heat. 

We  note  here  that  extinction  always  occurs  as  t  ^  oo.  However,  for  some  and  F, 
extinction  occurs  for  t  <  oo.  We  will  prove  the  former  in  Section  1.4. 

In  order  to  properly  consider  equation  (1.19)  so  that  we  may  determine  extinction 
times  for  various  and  F,  we  must  first  set  up  our  boundary  and  initial  conditions.  For 
our  boundary  condition,  we  will  choose  the  homogeneous  Dirichlet  condition, 

ii  =  0  on  X  [0,  oo).  (1-21) 

This  means  that  on  the  boundary  of  H,  dQ,  we  will  fix  the  temperature  to  zero  for  all  time 
t  >0.  We  could  have  fixed  the  boundary  to  be  some  nonzero  constant  temperature,  or  we 
may  also  have  chosen  the  temperature  at  the  boundary  to  change  over  time,  u  =  g{t),  for 
some  nonnegative  function  g. 
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Another  possible  choice  for  the  boundary  condition  is  a  Neumann  condition,  = 
g{t),  where  is  the  derivative  of  u  in  the  outward  normal  direction.  This  boundary 
condition  specifies  heat  flow  across  the  boundary  or  the  heat  flux.  We  could  have  also 
chosen  a  Robin  boundary  condition,  |^  +  Au  =  g{t),  in  which  there  is  an  interchange 
of  heat  between  the  edge  of  the  material  and  the  surrounding  medium.  Here  A  is  some 
constant. 

Now,  we  must  choose  our  initial  condition.  We  will  choose 

rt(x,0)  =  uo{x)  >0  on  H,  (1-22) 

where  H  H  U  dQ.  That  is,  uq  is  the  initial  temperature  of  the  material. 

Combining  equations  (1.19),  (1.21),  and  (1.22),  we  get  an  initial-boundary-value 
problem  (IB VP), 


Ut 

=  V2$(7x)  -  F{u) 

in  H  X  (0, oo). 

u 

=  0 

on  dQ  X  [0, oo), 

(1.23) 

u{x, 0) 

=  uoix)  >  0 

on  H. 

As  noted  earlier,  we  are  interested  in  the  extinction  time  for  solutions  to  equation 
(1.23).  However,  it  is  well-known  that,  in  general,  there  are  no  classical  solutions  to  this 
non-linear  parabolic  equation  for  arbitrary  choice  of  and  F.  In  Chapter  HI,  we  will  derive 
a  numerical  algorithm  to  calculate  the  solution  of  the  one-dimensional  version  of  equation 
(1.23).  In  Chapter  IV,  we  will  numerically  calculate  the  solution  for  several  choices  of  $ 
and  F  in  order  to  determine  whether  finite  extinction  time  occurs. 

Now,  let  us  look  at  a  related  problem  given  by  the  diffusion-free  equation, 

z'{t)  =  -[<^{z{t))  +  F{z{t))],  0<t<oo, 

z(0)  =  M  >  0. 
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It  can  be  shown  that  the  solution,  z{t),  to  equation  (1.24),  has  finite  extinction  time  if  and 
only  if 


ds 


/o  Hs)  +  F{s) 


<  oo. 


(1.25) 


where  e  >  0. 


We  want  to  determine  if  u  will  have  a  finite  extinction  time,  when  and  F  satisfy 
condition  (1.25).  So,  if  we  have  a  and  F  that  satisfy  condition  (1.25),  then  z  will  have 
a  finite  extinction  time.  If  we  then  compare  the  nonnegative  solutions  u  and  z,  and  find 
that  u{x,  t)  <  z{t)  for  all  x  and  t,  then  u  will  also  have  a  finite  extinction  time.  It  can  be 
shown  analytically  that  u  <  z  for  the  case  where  equation  (1.24)  has  no  term.  That  is, 
if  z  satisfies 


z'{t)  =  —F{z),  0  <t  <  oo, 

z(0)  =  M, 


(1.26) 


and  u{x,0)  <  M  where  tt  is  a  solution  of  (1.23),  then  u  <  z.  It  has  not  been  shown  for 
equation  (1.24). 


We  should  also  mention  that  if  (1.25)  holds,  and  we  can  calculate 


ds 


/o  Hs)+F{sy 


then  the  unique  solution  of  equation  (1.24)  is  z  where  z{t)  satisfies 

f 

JM 


rz{t) 

Im  ^{s)  +  F{s)  /o 


(1.27) 


(1.28) 


and  no  numerical  approximation  is  needed. 

The  numerical  results  for  equation  (1.24)  as  well  as  a  comparison  to  the  solution 
of  the  one-dimensional  version  of  (1.23)  are  given  in  Chapter  IV.  The  derivation  of  the 
numerical  solution  for  (1.24)  is  given  in  Chapter  III. 
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1.4  Proof  that  extinction  always  occurs  as  t  ^  oo 

Theorem  1.4.1  Suppose  u  is  a  nonnegative  solution  of  equation  (1.23),  where  ^  and  F 
are  nondecreasing,  nonnegative  (^^((OjOo))  nC'([0,oo))  functions  that  satisfy  (1.20).  Then 
In  F{s)dsdx  ^0  as  t  ^  oo. 

Proof.  First,  we  multiply  equation  (1.19)  by  F(u)  and  we  integrate  over  17  to  get 

/  F{u)utdx  =  /  F{u)V‘^^{u)dx  —  /  [F{u)]^dx.  (1-29) 

y  £7  </  £7  J  ^ 

Consider  the  term  J^F{u)V‘^^{u)dx.  We  see  that  with  integration  by  parts, 

[  F{u)V^^{u)dx  =  [  V[F(u)V^(u)]  -  [  VF{u)V^{u)dx.  (1.30) 

Jo.  Jq 

Now,  we  apply  the  divergence  theorem  to  equation  (1.30)  and  obtain 


in 


F{u)V‘^^{u)dx  =  /  F{u)V^{u)  ■  vds  —  /  VF{u)  •  V^{u)dx 

Jdn  Jn 


(1.31) 


=  —  /  VF{u)  ■  V^{u)dx, 

Jn 


where  n  is  the  unit  outward  normal  to  517,  ds  denotes  an  element  of  surface  area,  and 
F(u)V^(u)  ■  vds  =  0  because  u  vanishes  on  517  and  F(0)=0. 

Now,  substituting  (1.31)  back  into  equation  (1.29)  we  arrive  at 

/  F{u)utdx  =  —  VF{u)  ■  V^{u)dx  —  /  [F{u)]‘^dx 

J  Q,  J  ^  J  ^ 

=  -  [  F\u)<^\u)\Vu\‘^dx  -  [  [F{u)fdx  (1.32) 

Jn  Jn 

<  -  [  [F{u)fdx, 

Jn 

because  F{u)  and  <l>(u)  are  nondecreasing. 

Next,  we  define  G{u)  =  F{s)d.s  so  that  equation  (1.32)  becomes 

—  f  G{u{x,f))dx  <  —  f  [F{u{x,t))]'^dx.  (1.33) 

57  Jn  Jn 
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However,  since  F'{u)  >  0,  then 


ru  pu 

G{u)  =  /  F{s)ds  <  /  F{u)ds  =  uF{u). 
Jo  Jo 


(1.34) 


Taking  the  integral  of  both  sides  of  equation  (1.34),  we  get 


in 


G{u)dx  <  f  uF{u)  <  \  [  u'^dxi  f  [F{u)]‘^dx. 
Jn  \  Jn  \  Jn 


(1.35) 


On  a  side  note,  if  we  start  with  equation  (1.19),  and  multiply  both  sides  by  u,  then 
integrate  by  parts  and  apply  the  divergence  theorem  as  we  did  above,  we  obtain 


/  uutdx  =  /  uV‘^^{u)dx  —  /  uF(u)dx 

In  Jn  Jn 


=  —  Vu  •V^{u)dx  —  /  uF{u)dx 


=  —  Vu  •  ^'{u)Vu  —  /  uF{u) 


(1.36) 


=  -  $'(u)|Vup-  /  uF{u) 

Jn  Jn 

<  0, 


because  <h'(rt)  >  0,  F{u)  nonnegative,  and  u  >  0. 
We  see  that 


and  therefore. 


[  v?dx  =  2  [  uutdx, 
dt  Jn  Jn 


d  f 

—  /  u^dx  <  0. 

dtJn 


(1.37) 


(1.38) 


Equation  (1.38)  means  that  J^uJ{x,t)dx  is  bounded  above,  or  in  other  words. 


/  u‘^{x,t)dx<  /  UQ{x)dx  =  M^, 

In  Jn 


(1.39) 


where  uo{x)  is  the  initial  condition,  and  is  a  constant. 
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Now,  substituting  (1.39)  into  equation  (1.35)  we  obtain 


I  G{u)dx  <  mJ  j  [F{u)]^dx. 


(1.40) 


Squaring  both  sides  of  equation  (1.40)  and  dividing  by  — gives 


(1.41) 


Next,  we  combine  equations  (1.33)  and  (1.41)  to  get 


^  In  -~w(yl  G{u)dx'^  .  (1.42) 


Let  w{t)  =  J^G{u{x,t))dx,  so  that  the  inequality  in  (1.42)  becomes 


w\t)  <  -j^w\t). 


(1.43) 


Divide  both  sides  by  w‘^{t)  and  recognize  that 


w'{t)  _  d  (  I 


w‘^(t)  dt  \w{t) 


to  arrive  at 


dt  \w{t)  J 


(1.44) 


Next,  we  integrate  both  sides  with  respect  to  t  to  get 


1  1  t 

w{t)  t(;(0)  “  ' 


(1.45) 


Finally,  we  solve  for  w{t)  to  obtain 


0  <  w{t)  <  — j- 


_L  ^ 
0(0)  ^ 


(1.46) 


which  approaches  zero  as  t  ^  oo.  Thus,  extinction  always  occurs  as  t  ^  oo.  Q.E.D. 

In  summary,  we  know  that  all  nonnegative  solutions  of  (1.23)  have  extinction  time, 
though  it  may  be  at  infinity.  We  wish  to  determine  whether  extinction  occurs  in  finite 
or  infinite  time  for  any  given  and  F.  We  will  try  to  do  this  directly  by  numerically 
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calculating  the  solution  of  the  one-dimensional  version  of  equation  (1.23).  We  will  also 
try  to  do  this  indirectly  by  attempting  to  show  numerically  that  any  solution  u  of  the 
one-dimensional  version  of  (1.23)  must  satisfy  u  <  z  where  z  is  the  solution  to  (1.24)  with 
u{x,  0)  <  M. 
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II.  Literature  Review 


In  this  thesis,  we  are  studying  a  computational  model  of  non-linear  variables  and  F,  that 
predicts  temperatures  using  finite-difference  approximations.  In  predicting  temperatures, 
we  need  to  consider  specific  boundary  and  initial  conditions.  With  all  of  this  at  hand,  we 
hope  to  be  able  to  determine  whether  finite  extinction  time  occurs. 

In  this  chapter,  we  will  cite  a  few  examples  of  related  problems.  Among  these  are 
heat  pipe  problems,  dead  core  reaction-diffusion  problems,  and  non-linear  parabolic  finite 
extinction  time  problems. 

2.1  Heat  Pipe  Problems 

Heat  pipes  have  been  studied  for  many  years.  Among  these  studies  is  the  transient 
behavior  of  heat  pipes  under  many  external  conditions.  Under  normal  conditions,  it  has 
been  determined  that  transient  response  is  caused  by  the  thermal  capacity  and  conductance 
of  the  shell,  capillary  structure,  and  working  fluid,  and  is  only  slightly  influenced  by  liquid 
and  vapor  dynamics  [3]. 

Chang  and  Colwell  [3]  study  low-temperature,  heat  pipe  operation.  They  develop 
a  computational  model  for  predicting  the  transient  temperatures  of  low-temperature  heat 
pipes  based  on  finite-difference  approximations.  They  consider  variations  of  thermal  prop¬ 
erties  using  cubic  spline  interpolation  and  various  boundary  conditions  for  thermal  cou¬ 
pling.  They  also  calculate  nodal  temperatures  using  an  alternating  direction  implicit 
method. 

Chang  and  Colwell  then  compare  their  results  to  experimental  values  that  they  pro¬ 
duced  in  the  lab.  The  final  predicted  steady-state  values  are  in  good  agreement  with  the 
experimental  values.  However,  the  predicted  temperatures  reach  steady-state  faster  than 
the  measured  values  mainly  because  of  the  assumption  of  a  perfectly  insulated  heat  pipe 
and  uniform  vapor  temperature  [3]. 
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2.2  Dead  Core  Problems 


Dead  core  problems  are  very  interesting  in  nature.  The  problems  themselves  are 
brought  about  by  the  diffusion  of  one  substance  through  the  pores  of  a  solid  body  which, 
with  the  evolution  or  absorption  of  heat,  may  absorb  and  immobilize  some  of  the  diffusing 
substance.  The  heat  itself  will  then  diffuse  through  the  medium,  affecting  the  amount  of 
substance  the  solid  can  absorb  [4].  If  a  region  of  zero  reactant  concentration  is  formed  in 
finite  time,  we  have  a  dead  core  problem  [2].  More  specifically,  a  dead  core  is  a  region, 
Do  ^  D,  that  is  formed  in  a  finite  time  to  such  that  u{x,  t)  =  0  for  all  (x,  t)  G  Dq  x  [to,  oo). 

Bandle  and  Stakgold  [2]  study  non-linear  reaction-diffusion  problems  in  which  a  dead 
core  may  be  formed.  They  formulate  the  conditions  needed  for  a  dead  core  to  exist, 
and  calculate  estimates  for  its  time  of  onset.  The  simple  model  Bandle  and  Stakgold 
produce  is  a  system  of  homogeneous  non-linear  parabolic  equations  with  boundary  and 
initial  conditions.  The  main  parameters  of  which  are  both  concentration  and  temperature. 

Bandle  and  Stakgold  [2]  show  that  if  there  is  strong  reaction,  there  will  be  a  dead 
core  formation.  They  also  derive  estimates  for  the  time  of  onset,  size,  and  location  of  the 
dead  core.  Finally,  they  prove  convexity  of  the  dead  core. 

2.3  Finite  Extinetion  Time  Problems 

In  this  thesis,  we  are  concerned  with  the  non-linear  parabolic  finite  extinction  time 
problem  whose  solution,  u,  is  always  zero  on  the  boundary,  and  whose  initial  function,  uo, 
is  nonnegative.  Therefore,  we  will  consider  the  problem 


Ut 

=  A(p{u)  -  f{u) 

in  D  X  (0, oo), 

u 

=  0 

on  dD  X  [0,  oo 

u{x,  0) 

=  uo{x)  >  0 

on  D, 

where  D  is  a  bounded  domain  in  and  and  /  are  nondecreasing,  nonnegative 

(^^((O,©©))  n  (^([O,©©))  functions  that  satisfy 

(/?(0)  =  0,  /(O)  >  0,  (f{s)  >  0,  f{s)  >0  if  s  >  0.  (2.2) 
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Kalashnikov  [7]  and  Kersner  [8]  both  show  that  the  nonnegative  solution  of  the  one¬ 
dimensional  version  of  (2.1)  with  =  M,  ip"  >  0,  and  p'{^)  <  oo,  has  finite  extinction  time 
if  the  condition  of  strong  absorption, 


r#y<~.  (2.3) 

Jo  j{s) 

holds.  Conversely,  Kalashnikov  [7]  shows,  among  other  things,  that  if  (2.3)  is  infinite  and 
the  condition  p'(s)f{s)  <  Ks  for  some  constant  K  >  0  holds,  then  the  solution  does  not 
have  a  finite  extinction  time. 

Diaz  and  Diaz  [5]  consider  equation  (2.1)  with  /  =  0  and  show  that  if  D  is  bounded, 
then  a  necessary  and  sufficient  condition  for  u  to  have  finite  extinction  time  is  the  condition 
of  fast  diffusion,  i.e., 

<  oo.  (2.4) 

^(s) 

Lair  [9]  shows  that  if  either  (2.3)  or  (2.4)  holds,  then  the  solution  to  (2.1)  has  finite 
extinction  time.  Lair  and  Oxley  [10]  extend  Lair’s  result  to  show  that  if 

[  (  f(  \  =  (2-5) 

Jo  ns)  +  f{s) 

holds  then  there  is  no  finite  extinction  time.  Equivalently,  they  also  show  that  if  finite 
extinction  time  exists,  then  the  integral  in  (2.5)  is  finite. 

The  above  results  lead  to  an  open,  unsolved  problem.  Suppose  /  and  p  satisfy 


r  ds  _ 
Jo  f{s) 


ds 

p{s) 


=  oo, 


and 


ds 


V’is)  +  f{s) 


<  oo. 


(2.6) 


This  is  to  say  that  we  have  slow  diffusion  and  weak  absorption,  but  their  combined  effect  is 
strong.  Will  any  solution  of  equation  (2.1)  have  a  finite  extinction  time?  This  is  unknown 
even  for  the  special  case  where  p  and  /  are  both  concave  down  [10] . 

The  primary  objective  of  this  thesis  is  to  provide  an  answer  to  this  open  problem. 
Or,  if  unable  to  answer  it  directly,  to  provide  computational  evidence  of  the  correct  answer 
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and  perhaps  give  insight  into  how  to  prove  analytically  what  the  computations,  in  the 
examples  considered,  indicate. 
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III.  Numerical  Analysis 

In  this  chapter,  we  will  derive  our  numerical  method  for  the  solution  of  equation  (1.23)  in 
one-dimension.  We  will  use  a  predictor-corrector  method  in  order  to  ensure  the  accuracy 
of  our  solution.  Later  in  the  chapter,  we  will  derive  our  numerical  solution  to  equation 
(1.24). 


3.1  Corrector  for  absorption-diffusion  problem 


In  this  section,  we  will  develop  a  numerical  method  for  the  solution  of  the  one¬ 
dimensional  version  of  equation  (1.23)  with  H  =  (0,1).  We  start  by  choosing  a  method 
that  mirrors  the  Crank-Nicolson  method  because  when  this  method  is  used  to  solve  linear 
partial  differential  equations,  it  converges  and  is  numerically  stable  for  all  Ax  >  0  and 
At  >  0.  It  also  has  a  global  accuracy  of  O(Ax^)  [12].  Later,  we  will  see  that  we  need 
another  numerical  method  to  act  as  a  predictor  to  help  solve  this  method.  We  should  note 
here  that  this  implicit  method  can  be  made  into  a  system  of  linear  equations,  with  one 
equation  for  each  x  node.  This  system  can  then  be  written  in  a  matrix-vector  form  with 
a  tri-diagonal  matrix,  so  that  the  system  is  easily  solvable. 


Starting  with  equation  (1.18),  we  see  that  we  can  rewrite  this  in  terms  of  partial 
derivatives. 

ut  =  [4>(m)]xx  -  F{u) 


d  ,  .du 
dx 


-F{u) 


=  ^[^'{u)u.]-F{u) 

=  ^'{u)Ua;x  +  ^"{u){Ux)‘^ 


F{u). 


(3.1) 


Rewriting  the  equation  this  way,  allows  us  to  create  a  finite  difference  scheme.  We  could 
have  also  created  a  finite  difference  scheme  directly  from  equation  (1.18),  where 


*h('u)xa; 


‘h(u(x  -|-  Ax,  t))  —  24>(m(x,  t))  -|-  4>(ri(x  —  Ax,  t)) 

{A^ 


(3.2) 


but  this  method  was  found  to  be  unstable  for  4>(tt)  =  with  0  <  p  <  1. 
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First,  we  will  discretize  in  space.  Starting  with  the  Taylor  expansion  of  the  solution, 
u{x,t),  a  small  distance.  Ax,  away  from  x,  we  have 

A,x‘^  A 

u{x  +  Ax,  t)  =  u{x,  t)  +  Ax  Ux{x,  t)  H - —Uxx{x,  t)  H - —Uxxx{x,  t)-\ -  (3.3) 

2  0 


and 


u(x  -  Ax,  t)  =  u{x,  t)  -  Ax  Ux{x,  t)  H - —Uxx{x,  t) - —Uxxx{x,  t)-\ - ,  (3.4) 

2  b 

where  Ax”  =  (Ax)”.  Subtracting  (3.4)  from  (3.3)  and  truncating  the  series,  we  arrive  at 

A,x^ 

u{x  +  Ax,  t)  —  u{x  —  Ax,  t)  =  2Ax  Ux{x,  t)  H - —UxxxH,  t),  (3-5) 


where  x  —  Ax  <  ^  <  x  +  Ax. 

Rearranging  (3.5)  in  terms  of  Ux  we  get 


Ux{x,t) 


u{x  +  Ax,  t)  —  u{x  —  Ax,  t) 

2Ai 


Ax2 

~Y~ 


'^xxx 


(C,i). 


(3.6) 


Now,  in  our  particular  problem,  our  space  is  the  closed  interval  [0, 1] .  If  we  divide  this  space 
into  N  equal  intervals  of  length  Ax,  then  we  have  nodes,  Xj  =  iAx,  where  i  =  0, 1, . . . ,  A. 
Therefore,  xq  is  on  the  boundary  at  0  and  xm  is  at  1.  We  let  Ui{t)  =  u{xi,t)  in  equation 
(3.6)  and  we  regard  —^^^Uxxx{^:  t)  as  the  error  term  in  our  approximation  of  Ux-  We  note 
that  the  error  on  the  approximation  of  Ux  is  O(Ax^).  This  gives  us  an  approximation  for 
Ux  at  each  node. 


Ux{Xi,t) 


Ui+l{t)  -  Ui-l{t) 

2Ax 


(3.7) 


Going  back,  if  we  add  equations  (3.3)  and  (3.4),  truncate,  and  make  the  proper 
substitutions  as  described  above,  we  arrive  at  an  approximation  for  Uxx  at  each  node  Xj, 


'^xx{Xi,  t) 


Ui+i{t)  -  2ui{t)  + 
Ax^ 


(3.8) 


We  note  here  that  the  error  on  the  approximation  is  — ^Uxxxx{^,t)  which  is  O(Ax^). 
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If  we  substitute  our  newly  calculated  approximations,  (3.7)  and  (3.8),  into  equation 
(3.1),  we  arrive  at  a  spatially  discretized  version  of  our  problem. 


dt 


Ui{t) 


Ui+i{t)  -  2ui{t)  + 


+  ^>"(^^*(t)) 


Ui+l{t)  -  Ui-l{t)Y 

2Ax 


F{ui{t)). 


(3.9) 


At  this  point  we  see  that  we  have  derivatives  of  in  (3.9).  However,  since  we  know 
<h,  we  can  a  priori  calculate  its  derivatives,  and  substitute  them  in  where  necessary.  This 
seems  to  be  a  problem  for  cases  where  and  are  undefined  on  the  boundary  where 
=  0,  but  it  is  not  because  we  never  need  to  solve  for  a  solution  directly  on  any  boundary. 
This  is  because  our  boundary  conditions  force  our  solution  to  be  zero  on  the  boundaries. 

We  still  have  a  time  derivative  left  in  equation  (3.9),  and  therefore,  we  need  to 
discretize  in  time.  First,  we  will  take  the  integral  of  both  sides  of  equation  (3.9)  over  a 
small  interval  of  time  from  t  to  t  +  At. 


rt+At 


ut{xi,t)dt  = 


ft+At 


Ui+i{t)  -  2ui{t)  +  Ui-i{t) 


+  ‘h"(«,(t)) 


Ui+l{t)  -  Ui-l{t) 

2Ax 


-  F{ui{t))  >  dt. 


(3.10) 


Notice  that  ut{xi,  t)dt  =  Ui{t  +  At)  —  Ui{t).  Now,  with  the  remaining  integral  we  use 

the  trapezoid  rule  and  we  arrange  the  equation  so  that  only  Ui{t  +  At)  is  on  the  left  hand 
side  to  arrive  at 


Ui{t  +  At)  =  Ui{t)  +  ^  \  {ui{t  +  At)) 


Ui+i{t  +  At)  -  2ui{t  +  At)  +  Ui-i{t  +  At) 


+  ^"{ui{t  +  At)) 

+  d>'(rtj(t)) 

+  <^”{u^{t)) 


rtj+i(t  +  At)  -  Ui-i{t  +  At) 
2Ax 

Ui+i{t)  -  2ui{t)  +  Uj-iit) 

Ax^ 


Ax^ 

2 


-  F{ui{t  +  At)) 


Ui+l{t)  -  Ui-l{t) 

2Ax 


1  2 


-F{ui{t)) 


(3.11) 
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Here  we  stop  to  note  that  the  trapezoid  rule  on  a  linear  equation  has  a  local  inter¬ 
polation  error  of  [1].  This  means  that  our  use  of  the  trapezoid  rule  produces  a 

minimum  local  error  of  because  in  our  case  the  error  also  depends  on  F,  and 

our  O(Ax^)  approximations.  Also,  if  ‘h  =  tt,  |^|  is  small,  and  At  is  small,  then  this 
method  will  converge  [12].  It  is  our  assumption  that  this  method  will  converge  for  any 
or  F. 

In  our  problem,  the  time  interval  is  [0,  oo).  Dividing  this  interval  up  into  subintervals 
of  length  At,  we  arrive  at  nodes,  tn  =  nAt,  where  n  =  0, 1, . . .  .  This  means  to  is  at  the 
initial  time,  0.  We  will  now  let  Uin  =  Ui{tn)  =  u{xi,tn)  in  equation  (3.11)  thus  giving 


At 

+  ^"{Uin+l) 
+  ^'{Uin) 


J  (bUqi-  ,  ,i 

s  St'  yUifi+l) 

Ax^ 

'^i+ln+1  'W'j— In+l 

Ui-\-ln  2,Uin  +  Ui—ln 


A'('Ui  n+l) 


-h  ^"{Uin) 


2Ax 


-  F{uin)  I  . 

(3.12) 


Multiplying  out  the  squared  terms  and  reorganizing  the  terms  in  a  manner  similar 
to  the  arrangement  of  terms  in  the  Crank-Nicolson  method,  we  obtain 


Uin+l  1^1  T  ^2^2^  (^in+l) 

~^8Ax^*^  (^i  n-|-l)'W'j— 1  n+l 


+  Ui-^-ln+l 

T  Ui—l n+1 


At  / 

—  Uin 

+  ^^i+ln 

TgAx^^  ('*^*n+l)'U'i-|-l  n+1 


At  „  At  „ 

+  SAx^^  {uin)ui+ln  -  («in)wi-l) 

At  „  At  „ 

~*~8Ax2*^  ('*^*n)'Wi-ln  —  {Uin)ui+ii 


At  ,  At  „ 

2Ax2*^  SAx^^  (ifin+U^'-j+ln+l 

At  ,  At  „ 

'  2Ax^^  8Ax2^  ('Win+l)llj-ln+l 


+  Ui-li 


2Ax2 


{Uir 


2Ax2 


-  ^  [F{uin+l)  +  F{ui 


(3.13) 

Now  that  we  have  arranged  the  equation  the  way  we  want,  we  see  that  we  have  a 
numerical  scheme  for  our  problem.  However,  before  we  solve  our  problem,  we  can  make 
our  work  a  little  easier  by  writing  equation  (3.13)  in  matrix  form.  This  will  allow  us  to 
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solve  the  entire  x  interval  at  once  for  each  time  step. 


1  0  0 

Cln+l  O'ln+1  biri+1 
0  C2n+1  0.2  n+1 

0  0  0 

0  0  0 

1  0  0 

fin  din  Oin 

0  f2n  d2n 

0  0  0 

0  0  0 


0  0  'UOn+1 

0  0  rtin+1 

0  0  U2n+1 

ON—ln+1  bj\f—in+l  "ON—ln+l 

0  1  UNn+1 


0 

0 

0 

dN-ln 

0 


F{uon+l)  +  F{uon) 

F(^Uln+l)  F  F{u\n^ 

At  F{u2n+l)  +  F{u2n) 

~2  : 

F{uN-ln+l)  +  F{UN-In) 
F{uNn+l)  +  F{uNn) 


(3.14) 
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where 


n+l 
bi  n+l 
Q  n+l 
di  n 

n 

fi  n 


^  Ax2^  ('W'jn+l)) 

At  ,  At  n  At  n 

~  2  Ax^^  “  SAx^^  \Ain+l)Ui+ln+l  +  ^  y^in+ljUi-ln+l-, 

At  ,  At  ,  At  , 

~2Ax2^  ('*^*^+1)  “  SAx^^  \Uin+l)Ui-ln+l  +  (Ui  n+l)Ui+l  n+l  ■, 

^  ~  A^*^  {'^in), 

At  ,  At  ,  At  , 

2Ax2^  +  8Ax2^  “  8Ax2^  («in)Mi_ln, 

2Ax2^  +  8Ax2^  “  8Ax2*^  («in)Mi+l„. 


Note  that  we  can  incorporate  our  boundary  conditions  of  won  =  0  and  UNn  =  ^  into 
equation  (3.14).  This  will  simplify  the  equation  to 


ai  n+l 

bl  n+l 

0 

0 

111  n+l 

C2n+1 

0.2  n+l 

fe2n+l 

0 

ll2n+l 

0 

C3  n+l 

as  n+l 

0 

ll3n+l 

0 

0 

0 

UTV— ln+1 

llAf-ln+l 

dl  n 

^1  n 

0  • 

0 

111  n 

F{ui  n+l. 

+  F{uin) 

f2n 

d2  n 

62  n 

0 

112  n 

At 
~  ~2 

F{u2  n+l. 

+  F{u2n) 

0 

fsn 

CO 

0 

113  n 

F{u;^  n+l. 

-h  F{u3n) 

0 

0 

0  • 

1 

s 

UN-ln 

F{uN-ln+l) 

+  F{uj,f-in) 

(3.15) 

where  ain+i-,  hin+i-,  Cin+i-,  din,  Gin,  and  fin  are  the  same  as  in  equation  (3.14).  We  should 
stop  to  note  that  this  method  has  a  minimum  global  error  of  O(Ax^)  because  it  degenerates 
into  the  linear  Crank-Nicolson  method  as  will  be  shown  in  Section  3.2.  The  exact  error  is 
dependent  upon  the  choice  of  4>  and  F. 
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Representing  the  ain+i,  hn+i,  Qn+i  matrix  in  equation  (3.15)  by  ^n+i;  the  Uin+i 
vector  with  Un+i,  the  din,  ^in,  fin  matrix  with  Bn,  the  Uin  vector  with  Un,  and  the 
F{uin+i)  +  F{uin)  vector  with  Fn,  we  obtain 

■^n+lUn+l  —  Fntln  ~l^^n-  (3.16) 

Finally,  we  solve  for  Un+i  to  get 

Un+l  =  ^n+l^'n-^n  ^n+l^n-  (3-17) 

When  solving  equation  (3.17),  we  run  into  a  problem.  Unlike  the  linear  case  of  the 
Crank-Nicolson  method,  we  need  to  estimate  the  solution  at  the  n  + 1  time  step  in  order  to 
calculate  the  solution  at  the  n  +  1  time  step.  This  occurs  because  the  non-linearity  of  our 
equation  adds  extra  Uin  and  Uin+i  terms  to  the  tri-diagonal  matrices.  However,  this  is  not 
a  major  problem  since  we  can  use  equation  (3.17)  as  our  corrector  in  a  predictor-corrector 
method. 

We  will  choose  a  forward  difference  method  as  our  predictor  because  it  is  a  good, 
standard  explicit  method  with  a  local  error  of  O(Ax^)  for  linear  problems  [1].  By  making 
a  good  guess  at  our  solution  with  the  forward  difference  method,  we  have  assured  that 
our  solution  has  a  minimum  global  error  of  O(Ax^).  This  is  particularly  important  when 
solving  a  non-linear  problem  such  as  ours.  The  predictor  will  be  formulated  in  Section  3.3. 

3.2  Check  of  corrector  formulation 

Before  we  move  on  to  the  formulation  of  the  predictor,  we  should  first  check  to  see  if 
our  method  degenerates  to  the  Crank-Nicolson  method  when  our  problem  becomes  linear. 
We  must  realize  that  for  the  linear,  one-dimensional  case  of  equation  (1.23),  ‘h(tt)  =  u. 
This  means  that  ^'{u)  =  1  and  =  0.  Also,  the  absorption  term  will  only  depend 

upon  X  and  t  in  the  linear  case.  Thus,  F{u{x,  t))  =  F(x,  t).  With  all  this  in  mind,  we  shall 
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substitute  these  values  into  (3.15)  to  give  us 


’  1  +  ^ 

^  ^  Aa? 

-At 

2Ax2 

0 

0 

Uln+l 

-At 

1  +  ^ 

^  ^  Ax2 

-At 

0 

U2n+1 

2Ax'^ 

2Ax2 

0 

-At 

2Ax2 

^  ^  Ax2 

0 

U3n+1 

0 

0 

0 

^  Ax2  J 

UN-ln-el 

I  At 

Ax^ 

At 

2Ax2 

0 

0 

" 

Ul  n 

At 

1  At 

Ax2 

At 

0 

U2n 

2Ax2 

2Ax2 

0 

At 

2Ax2 

At 

Ax2 

0 

U3n 

0 

0 

0 

1 

UN-ln 

At 


F(^Xi,  -|-  F(^Xi,  tn) 

F{x2,tn+l)  +  F{x2,tn) 
F{x^,  -|-  i^(x3,  tfi) 

F{xN-l,tn+l)  +  F{XN-I,tn) 


(3.18) 


which  is  the  Crank-Nicolson  method  for  a  linear  absorption-diffusion  equation  with  zero 
Dirichlet  boundary  conditions. 

We  note  here  that  by  degenerating  the  corrector  method  to  the  Crank-Nicolson 
method,  and  if  F  =  0,  we  remove  the  need  for  a  predictor.  This  is  because  there  are 
no  longer  any  Uin+i  terms  in  the  matrices  of  equation  (3.18).  The  method,  however,  is 
still  an  implicit  method. 


3.3  Predictor  for  absorption-diffusion  problem 

We  have  already  decided  to  use  a  forward  difference  method  as  our  predictor  so  we 
begin  as  before  with  equation  (1.18).  We  then  follow  all  of  the  same  steps  as  in  Section  3.1 
until  we  arrive  at  equation  (3.10).  Again  we  note  that  ut{xi,  f)dt  =  Ui{t  +  At)  —  Ui{t), 
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and  we  make  the  substitution  to  give 


Ui{t  +  At)  -  Ui{t) 


rt+At 


'ui+i{t)  -  2ui{t)  +  Ui-i{ty 

Ax^ 

Ui+l{t)  -  Ui-l{t) 

2Ax 


-  F{ui{t))  >  dt. 


(3.19) 


We  will  now  use  the  forward  rectangle  rule  on  the  remaining  integral  in  equation 
(3.19).  We  will  also  simultaneously  make  the  substitution  Uin  =  Ui{tn)  and  leave  only  the 
Uin+i  term  on  the  left  hand  side.  This  gives 


Uin+l 


Ui 


+  At  |^>'(iiin) 


Ui-\-ln  ‘2‘Uin  T 
Ax^ 


+  {Uin) 


Ui—\n 

2Ax 


F(^Ui  1 


(3.20) 


It  is  from  this  step  that  we  gain  a  truncation  error  which  has  a  minimum  error  of  0{At^). 
The  actual  error  depends  on  our  choice  of  and  F. 

Multiplying  out  all  the  terms  in  (3.20)  and  rearranging  in  an  appropriate  manner, 
we  arrive  at 


n+1 


1  2At  ■ 


T  ^i+ln 

T  Ui—ln 


At  ,  At  , 

^  (^in)  T  Tyr  2^  (^in)^i+l) 


Ax^ 


dAx^ 
At 


-  AtF{uin)- 


(3.21) 

Now,  we  have  a  numerical  scheme  for  our  problem,  but  just  as  in  Section  3.1,  we  can  make 
our  work  a  little  easier  by  writing  equation  (3.21)  in  matrix  form.  This  will  allow  us  to 


3-9 


solve  the  entire  x  interval  at  once  for  each  time  step. 


UOn+l 

1 

0 

0  • 

0 

0 

"U-On 

Uln+1 

Pin 

9l  n 

hi  n 

0 

0 

'W'l  n 

U2  n+1 

= 

0 

P2n 

92  n 

0 

0 

U2n 

UN-ln+1 

0 

0 

0 

9N-ln 

hN-ln 

UN—1  n 

UNn+l 

0 

0 

0  • 

0 

1 

UNn 

-  At 


F{uon) 

F{uin) 

F{u2n) 

F{UN-In) 

F{uNn) 


(3.22) 


where 


9in  — 
hin  — 
Pin  — 


1  .  r,  {Ui  n) : 


At 

Ax^ 

At 


At 


^  i^in)  F  ^^^2^  n)'U'i+l  n 


At 

4Ax^ 


^  (Uin)Ui—lnj 


4Ax^^  (^^in)Wi-ln  4Ax^^  (Uin)Ui+lt 


We  can  now  incorporate  our  boundary  conditions,  rton  =  0  and  ujvn  =  0,  into 
equation  (3.22).  This  will  simplify  the  equation  to 


^tln+l 

9l  n 

hi  n 

0  • 

0 

'^1  n 

F{uin) 

U2n+1 

P2n 

92  n 

h2n 

0 

^2  n 

F{u2n) 

UZn+l 

= 

0 

P3  n 

CO 

0 

'^3  n 

-  At 

F{u3n) 

UN—1  n+1 

0 

0 

0  • 

1 

UN-ln 

F{UN-In) 

(3.23) 
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where  gim  him  and  pin  are  the  same  as  in  equation  (3.22). 

Representing  the  Uin+i  vector  in  equation  (3.23)  with  the  gim  him  Pin  matrix 

with  Cm  the  Uin  vector  with  Un  \  and  the  F{uin)  vector  with  F,  we  obtain 

-  AtP,  (3.24) 

where  the  superscript  (p)  denotes  predictor. 

We  now  have  enough  information  to  generate  an  accurate  numerical  solution  to  our 
problem.  Before  we  do,  we  should  note  that  equation  (3.24)  has  the  stability  condition 
that  <  5-  We  shall  show  this  in  the  next  section. 

3-4  Stability  condition  for  the  predictor  of  the  absorption- diffusion  problem 

In  order  to  show  that  our  predictor  has  the  stability  condition  ^ <  5,  we  begin 
with  equation  (3.1), 

ut  =  ^'{u)uxx  +  -  F{u).  (3.25) 

In  general,  the  absorption  term  plays  no  role  in  stability,  so  we  will  neglect  this  term.  Now, 
we  make  a  first  order  approximation  for  ^”{u){uxp,  which  is  zero.  This  leaves  us  with 

Ut  =  ^\u)uxx,  (3.26) 

which  contains  the  major  contributors  to  the  predictor  methods  stability  condition. 

Discretizing  equation  (3.26)  in  the  same  manner  as  the  predictor  method  described 
in  Section  3.3  gives 

^in+1  —  ^in  /\x‘^  ‘  (3.27) 

We  will  now  use  von  Neumann’s  method.  First,  we  get 

pn+i^ij9  ^  ^  -  2p^F^^  +  ,  (3.28) 

L  J 
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where  j  =  \/— 1  and  6  is  arbitrary. 


We  solve  for  our  error  growth  factor,  p,  to  obtain 

^\u)At 


p  =  1  + 


=  1  +  2 


Ax^ 


ejo  +  -  2 

[cos(0)  —  1] 


Ax2  V2 


Now,  we  need  |p|  <  1.  Therefore, 


At  2  f  0\ 

-|<1, 


which  with  some  manipulation  becomes 


Since  0  is  arbitrary. 


where  <!>'  >  0. 


0  < 


^'{u)At 

Ax^ 


sin^ 


<!>'{u)At  ^  1 
Ax2  -  2’ 


(3.29) 


(3.30) 


(3.31) 


(3.32) 


3.5  Use  of  predictor-corrector  for  absorption- diffusion  problem 

Now,  we  just  need  to  put  it  all  together.  First,  we  choose  a  Ax  and  At  that  satisfy 
the  stability  condition,  ^ <  5,  of  the  predictor.  By  doing  this,  we  assure  that  our 
predictor  method  will  be  computationally  stable.  This  means  that  if  we  introduce  an  error 
at  any  time  step,  the  error  will  not  increase  as  we  proceed.  Next,  we  know  that  our  initial 
condition  is  tt(x,0)  =  uo(x),  so  we  create  our  initial  vector,  from  this.  We  then  solve 
equation  (3.24)  for  Now,  we  substitute  and  into  Bn,  Un,  and  Fn  in 

equation  (3.16).  Next,  we  solve  (3.16)  for  ui  as  shown  in  equation  (3.17).  Then,  we  let 
ui  =  and  we  use  this  as  our  new  predictor  in  equation  (3.24).  Finally,  we  repeat  the 
process  for  each  time  step  until  we  have  reached  our  final  time  step.  In  cases  with  finite 
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extinction  time,  our  final  time  step  is  at  the  time  of  extinction.  For  cases  with  infinite 
extinction  time,  we  allow  the  computer  to  run  until  Uin  <  See  Section  4.2  for  a 

more  detailed  explanation. 


3.6  Predictor-corrector  formulation  and  use  for  diffusion-free  problem 

Now  that  we  have  produced  a  numerical  method  which  can  solve  the  one-dimensional 
version  of  equation  (1.23)  with  zero  Dirichlet  boundary  conditions,  we  can  move  on  to 
creating  a  numerical  method  for  solving  equation  (1.24).  We  shall  start  by  taking  the 
integral  of  both  sides  of  equation  (1.24)  from  time  t  to  t  At  in  order  to  get 

Jt 


/t-cAt 

[4>(z(t))  -|-  F{z{t))]dt. 


(3.33) 


We  note  that  z'{f)dt  =  z{t-\-  At)  —  z{f)  and  so  we  substitute  this  into  equation  (3.33) 


in  order  to  get 


/t+At 

[^{z{t))  -I-  F{z{t))]dt. 


(3.34) 


We  now  use  the  trapezoid  rule  to  evaluate  the  integral  on  the  right  hand  side  of 
(3.34).  This  gives 


z{t  -I-  At)  -  z{t)  =  -^[^{z{t  -I-  At)  -I-  F{z{t  -I-  At))  +  4>(2;(t))  -|-  F{z{t))].  (3.35) 


The  truncation  error  from  using  the  trapezoid  rule  is  at  minimum  0{At^)  and  depends  on 
the  choice  of  and  F. 

Now,  we  can  divide  our  time  interval,  [0,  oo),  into  subintervals  of  length  At  to  obtain 
nodes,  tn  =  nAt,  where  n  =  0, 1, . . .  .  This  means  to  =  0.  It  also  means  that  we  can 
make  our  nodes  the  same  distance  apart  as  we  have  done  for  our  numerical  solution  of  the 
one-dimensional  version  of  (1.23)  just  by  choosing  the  same  At.  This  will  allow  us  to  more 
easily  compare  the  solutions  of  equations  (1.23)  and  (1.24).  We  let  Zn  =  z{tn)  in  equation 
(3.35)  and  we  separate  the  n  and  n  -|-  1  terms  thus  giving 

Zn+l  +  ^[^{Zn+l)  +  F{Zn+l)]  =  Zn  +  ^[^{Zn)  +  F{Zn)]-  (3.36) 
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Now,  we  have  an  implicit  method  for  the  solution  of  (1.24).  This  means  that  we 
will  have  to  make  an  estimate  for  the  solution  at  the  n  +  1  time  step  in  order  to  solve 
equation  (3.36).  We  will  do  this  with  a  predictor,  and  we  will  make  (3.36)  our  corrector. 
For  our  predictor,  we  will  go  back  to  equation  (3.34).  Instead  of  using  a  trapezoid  method 
to  estimate  the  integral  on  the  right  hand  side  of  the  equation,  we  will  use  a  forward 
difference  method.  By  doing  this  we  get 

z{t  +  At)  —  z{t)  =  —  At[4>(2;(t))  +  F(2;(t))].  (3.37) 

This  step  adds  a  minimum  truncation  error  of  O(At^)  depending  on  $  and  F. 

Again,  we  make  our  Zn  =  z{tn)  substitution  and  solve  for  Zn+i  to  give  us 

Zn+l  =  Zn-  At[^{Zn)  +  F{Zn)].  (3.38) 

As  in  Section  3.3,  we  now  change  the  notation  slightly  so  that  we  can  follow  the  same 
procedure  as  in  Section  3.5  to  numerically  solve  equation  (1.24)  with  a  predictor-corrector 
method.  So,  equation  (3.38)  is  now 

=  4'”  -  A([3>(4''>)  +  (3.39) 

It  occurs  to  us  at  this  point,  we  could  have  just  used  the  predictor  method  as  our  sole 
solution  to  equation  (1.24).  Using  just  the  predictor  would  give  us  a  minimum  local  error 
of  0{At)  which  is  accurate,  but  when  dealing  with  a  non-linear  problem  such  as  ours,  we 
wish  to  be  as  accurate  as  possible.  So,  we  use  the  predictor-corrector  method  to  increase 
our  accuracy. 
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IV.  Numerical  Results 


In  this  chapter,  we  will  see  that  our  numerical  method  for  the  solution  of  the  one-dimen¬ 
sional  version  of  equation  (1.23)  is  an  accurate  solution  when  the  problem  is  linear.  We 
will  also  see  that  our  method  can  show  that  for  certain  cases  of  and  F,  the  solution  has 
finite  extinction  time,  consistent  with  the  analytical  results  of  Lair  [9].  At  the  end  of  the 
chapter,  we  will  compare  the  solutions  of  equation  (1.23)  and  equation  (1.24). 


4-1  Accuracy  of  the  numerical  solution  of  the  absorption-diffusion  problem 

We  know  that  our  numerical  method  for  the  one-dimensional  version  of  equation 
(1.23)  degenerates  to  the  Crank-Nicolson  method  for  linear  problems.  We  also  know  that 
for  linear  problems,  the  Crank-Nicolson  method  has  an  O(Ax^)  global  accuracy.  This 
means  that  we  can  test  a  linear  problem  with  a  known  analytic  solution  and  compare  it 
with  our  numerical  results. 

First,  we  will  choose  an  equation  with  a  known  solution.  Let  us  start  with  a  simple 
diffusion  equation  problem  with  zero  boundary  conditions  and  an  initial  function  uo{x)  = 
sin(7rx).  So,  we  shall  solve 


ut{x,f) 
u(0,  t) 
u{l,f) 
u{x,  0) 


Uxx{x,f),  0  <  X  <  1, 

0, 

0, 

sin(7rx),  0  <  x  <  1. 


0  <  t  <  oo, 
0  <  f  <  oo, 
0  <  f  <  oo. 


(4.1) 


We  know  that  this  has  an  exact  solution  of  Ue{x,t)  =  u{x,t)  =  e  sin(7rx). 

Now  that  we  have  chosen  a  problem  with  an  exact  solution,  we  will  now  choose  Ax 
and  At.  We  will  start  by  choosing  multiple  values  of  Ax.  This  way,  we  will  be  able  to 
compare  each  solution  to  the  exact  solution.  This  will  allow  us  to  see  if  we  achieve  the 
correct  order  of  accuracy.  We  will  choose  Ax  equal  to  0.1,  0.05,  and  0.02. 

Even  though  we  do  not  need  to  satisfy  the  stability  requirement  for  the  predictor 
because  this  a  linear  problem  with  F  =  0  (see  Section  3.2),  we  will  anyway.  So,  we  must 
choose  At  to  satisfy  the  condition  <  I  because  4>'(m)  =  1.  Therefore,  At  must  be  less 
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than  or  equal  to  0.005,  0.00125,  and  0.0002,  respectively.  We  will  choose  At  =  0.001  for 
Ax  =  0.1  and  0.05,  and  we  will  choose  At  =  0.0002  for  Ax  =  0.02. 

We  will  make  calculations  at  t  =  0.5,  1,  and  5.  For  ease  of  use  and  lack  of  space,  we 
only  show  the  solutions  at  spatial  increments  of  Ax  =  0.1.  We  will  define  our  solutions  for 
Ax  equal  to  0.1,  0.05,  and  0.02  to  be  u,i{x,t),  u,05{x,t),  and  u,02{x,t),  respectively.  Also, 
we  will  define  our  errors  to  be  e,i(x,  t)  =  u,i{x,  t)  —  Ue{x,  t),  e.05(x,  t)  =  u,o^{x,  t)  —  Ue{x,  t), 
and  em{xjt)  =u,02{x,t)  —Ue{x,t). 


Table  4.1.  Numeric  and  exact  results  of  equation  (4.1)  at  t  =  0.5.  For  Ax  =  0.1  and 
0.5,  At  =  0.001.  For  Ax  =  0.02,  At  =  0.0002.  All  values  are  multiplied  by 
10-3. 


u(.l,.5) 

2.314072 

n(.2,.5) 

4.401627 

u(.3,.5) 

6.058320 

u(.4,.5) 

7.121982 

u(.5,  .5) 

7.488495 

u(.6,  .5) 

7.121982 

u(.7,.5) 

6.058320 

bo 

4.401627 

u(.9,  .5) 

2.314072 

U.02 

2.244971 

2.226021 

4.270189 

4.234144 

5.877411 

5.827799 

6.909311 

6.850989 

7.264879 

7.203556 

6.909311 

6.850989 

5.877411 

5.827799 

4.270189 

4.234144 

2.244971 

2.226021 

Ue 

e.i 

2.222414 

0.091658 

4.227283 

0.174344 

5.818356 

0.239964 

6.839888 

0.282094 

7.191883 

0.296612 

6.839888 

0.282094 

5.818356 

0.239964 

4.227283 

0.174344 

2.222414 

0.091658 

e.05 

e.02 

0.022557 

0.003607 

0.042906 

0.006861 

0.059055 

0.009443 

0.069423 

0.011101 

0.072996 

0.011673 

0.069423 

0.011101 

0.059055 

0.009443 

0.042906 

0.006861 

0.022557 

0.003607 

We  see  from  Table  4.1,  that  at  t  =  0.5,  we  improve  our  accuracy  as  we  lower  our 
Ax,  just  as  we  should.  If  we  average  the  errors  for  each  node  of  a  given  Ax,  we  get  the 
average  error  for  the  numerical  solution  with  that  Ax.  With  Ax  =  0.1,  our  average  error 
is  2.08081  X  10“^.  For  Ax  =  0.05  we  get  an  average  error  of  5.12087  x  lO-^,  and  we  get 
an  average  error  of  8.18856  x  10“®  for  Ax  =  0.02. 

In  going  from  Ax  =  0.1  to  Ax  =  0.05,  we  double  the  number  of  nodes.  If  our  method 
is  actually  O(Ax^),  then  our  error  should  be  reduced  by  a  factor  of  4.  If  we  divide  our 
average  error  at  Ax  =  0.1  by  4,  we  get  5.20203  x  lO-^,  which  is  approximately  equal  to  our 
average  error  with  Ax  =  0.05.  Therefore,  we  see  that  our  numerical  method  for  equation 
(4.1)  has  a  global  error  of  O(Ax^). 

We  can  also  look  at  the  relative  error.  Relative  error  is  defined  as  — ,  where  i  =  0.1, 
0.05,  and  0.02.  If  we  average  the  relative  errors  of  the  nodes  for  each  Ax,  we  get  the 
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average  relative  error  for  that  Ax.  Our  average  relative  errors  at  t  =  0.5  are  0.04124, 
0.01015,  and  0.00162,  for  Ax  =  0.1,  0.05,  and  0.02,  respectively. 


Table  4.2.  Numeric  and  exact  results  of  equation  (4.1)  at  t  =  1.  For  Ax  =  0.1  and  0.5, 
At  =  0.001.  For  Ax  =  0.02,  At  =  0.0002.  All  values  are  multiplied  by  10“®. 


u.i 

^^.05 

U.02 

Ue 

e.i 

e.05 

e.02 

ii(.l,l) 

1.732892 

1.630945 

1.603527 

1.598334 

0.134558 

0.032611 

0.005193 

u(.2,l) 

3.296156 

3.102241 

3.050089 

3.040213 

0.255943 

0.062028 

0.009876 

«(.3,1) 

4.536770 

4.269868 

4.198087 

4.184494 

0.352276 

0.085374 

0.013593 

u{A,l) 

5.333292 

5.019531 

4.935148 

4.919167 

0.414125 

0.100364 

0.015981 

u{.5,l) 

5.607756 

5.277847 

5.189121 

5.172319 

0.435437 

0.105528 

0.016802 

ii(.6,l) 

5.333292 

5.019531 

4.935148 

4.919167 

0.414125 

0.100364 

0.015981 

u{.7,l) 

4.536770 

4.269868 

4.198087 

4.184494 

0.352276 

0.085374 

0.013593 

u{.8,l) 

3.296156 

3.102241 

3.050089 

3.040213 

0.255943 

0.062028 

0.009876 

u{.9,l) 

1.732892 

1.630945 

1.603527 

1.598334 

0.134558 

0.032611 

0.005193 

Calculating  the  average  errors  at  t  =  1,  we  arrive  at  3.05471  x  10“®,  7.40313  x  10~^, 
and  1.17876  x  10“'^  for  Ax  =  0.1,  0.05,  and  0.02,  respectively.  Again  we  see  that  our 
global  error  is  of  O(Ax^).  Now  we  calculate  the  average  relative  errors  to  obtain  0.08419, 
0.02040,  and  0.00325,  respectively. 


Table  4.3.  Numeric  and  exact  results  of  equation  (4.1)  at  t  =  5.  For  Ax  =  0.1  and  0.5, 
At  =  0.001.  For  Ax  =  0.02,  At  =  0.0002.  All  values  are  multiplied  by  10“^^. 


u.i 

um 

u.02 

Ue 

e.i 

e.05 

e.02 

ii(.l,5) 

1.713672 

1.265513 

1.162657 

1.143954 

0.569718 

0.121559 

0.018703 

u(.2,5) 

3.259598 

2.407148 

2.211505 

2.175931 

1.083667 

0.231217 

0.035574 

u(.3,5) 

4.486452 

3.313155 

3.043875 

2.994912 

1.491540 

0.318243 

0.048963 

«(.4,5) 

5.274141 

3.894847 

3.578290 

3.520730 

1.753411 

0.374117 

0.057560 

u{.5, 5) 

5.545560 

4.095285 

3.762436 

3.701914 

1.843646 

0.393371 

0.060522 

u{.6, 5) 

5.274141 

3.894847 

3.578290 

3.520730 

1.753411 

0.374117 

0.057560 

u{.7,5) 

4.486452 

3.313155 

3.043875 

2.994912 

1.491540 

0.318243 

0.048963 

u{.8, 5) 

3.259598 

2.407148 

2.211505 

2.175931 

1.083667 

0.231217 

0.035574 

u{.9, 5) 

1.713672 

1.265513 

1.162657 

1.143954 

0.569718 

0.121559 

0.018703 

The  average  errors  at  t  =  5  are  1.29337  x  10“^^,  2.75960  x  10“^^,  and  4.24580  x  10“^^, 
respectively.  The  average  relative  errors  are  0.49802,  0.10626,  and  0.01635,  respectively. 
We  notice  that  our  relative  errors  seem  to  grow  in  an  almost  linear  relation  to  t.  This 
means  our  solution  will  eventually  become  relatively  inaccurate.  At  t  =  5  with  Ax  =  0.1, 
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we  are  already  starting  to  see  a  loss  of  relative  accuracy.  So,  it  is  important  to  choose  Ax 
and  At  wisely. 

Since  this  is  the  simplest  one-dimensional,  linear  case  of  the  diffusion  equation,  our 
numerical  method  will  generate  the  minimum  amount  of  error  possible.  So,  thinking  ahead, 
we  may  believe  we  could  have  accuracy  problems  for  large  t,  especially  for  poorly  chosen 
Ax  and  At.  Since,  in  general,  there  are  no  classical  solutions  to  (1.23),  we  cannot  test  the 
accuracy  of  our  method  for  non-linear  equations.  The  only  test  that  we  can  do  is  to  try 
various  and  F  to  see  if  finite  extinction  time  occurs,  and  that  these  results  are  consistent 
with  those  of  Lair  [9].  This  is  done  in  the  next  section. 

4-2  Finite  extinction  time  of  solutions 

Now  that  we  know  our  numerical  solution  is  accurate,  at  least  for  linear  equations, 
we  need  to  see  if  it  can  correctly  predict  finite  extinction  time.  To  do  this,  we  will  use  a 
modification  of  Lair’s  [9]  result.  We  consider  the  equation 

ut  =  [u^jxx  -  (4.2) 

for  some  p  >  0  and  q  >  0.  We  know  that  the  solution  of  (4.2)  has  finite  extinction  time 
if  and  only  if  min{p,  g}  <  1.  For  all  of  the  calculations  in  this  section  we  assume  the  zero 
Dirichlet  boundary  condition  and  the  initial  condition, 

ii(0,  t)  =  0,  0  <  t  <  oo, 

u{l,t)  =  0,  0  <  t  <  oo, 

tt(x, 0)  =  sin(7rx),  0  <  x  <  1. 

First,  we  will  test  the  case  where  p  =  ^  and  q  =  1.  This  corresponds  to  =  ^/u  and 

F  =  u.  This  case  should  have  finite  extinction  time,  and  we  see  from  Figure  4.1  that  the 
solution  has  an  extinction  time  at  t  =  0.205. 

This  is  a  good  time  to  stop  and  note  that  the  code  used  to  calculate  the  above 
solution,  and  all  of  the  other  solutions  in  this  chapter,  has  a  command  built-in  that  will 
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Figure  4.1.  Solution  to  equation  (1.23)  with 
At  =  0.001. 


1 


not  allow  the  calculated  solution  to  become  negative.  If  the  calculated  solution  at  any 
node,  for  any  iteration,  becomes  negative,  the  code  will  set  the  calculated  solution  at  that 
node  to  zero.  Therefore,  once  the  calculated  solution  at  any  node  reaches  zero,  it  cannot 
go  any  lower.  It  can,  however,  go  back  up  again  from  zero,  but  this  was  never  observed. 
The  justification  for  this  command  is  that  we  are  only  interested  in  nonnegative  solutions 
to  our  equation. 

We  shall  also  note  that  this  is  a  case  where  ‘h'  and  are  undefined  at  x  =  0,  in 
fact,  lim^^o =  co  and  lims^o =  oo.  We  choose  Ax  =  0.1  and  At  =  0.001 
because  they  satisfy  the  stability  condition  <  2  at  t=0,  and  the  computer  is  able 

to  handle  the  amount  of  processing  needed  to  calculate  the  numerical  solution  under  these 
conditions.  With  Ax  =  0.05,  we  need  to  choose  a  At  so  small  that  the  computer  cannot 
handle  the  processing  needed  to  calculate  the  numerical  solution. 
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For  the  next  few  examples,  we  will  choose  different  combinations  of  and  F  so  that 
our  solutions  will  have  finite  extinction  times.  We  shall  choose  =  y/u  and  F  =  for 
our  second  example.  This  particular  solution  reaches  extinction  at  t  =  0.212  as  shown  in 
Figure  4.2.  Again  we  note  here  that  we  have  chosen  Ax  =  0.1  and  At  =  0.001  because  of 
the  computing  problems  with  smaller  Ax  increments. 


0.25  0 

Figure  4.2.  Solution  to  equation  (1.23)  with  $  =  y/u  and  F  =  v^.  Ax  =  0.1  and 
At  =  0.001. 


We  see  that  by  making  q  bigger  and  keeping  p  the  same,  we  cause  extinction  to 
occur  at  a  slightly  later  time,  but  only  very  slightly.  In  our  next  two  examples,  we  will 
look  at  what  happens  when  we  keep  q  constant,  and  increase  p.  We  will  choose  =  n  and 
F  =  y/u,  and  and  F  =  y/u.  Figure  4.3  shows  that  extinction  occurs  at  t  =  0.468  for 

the  first  case,  and  for  the  second  case.  Figure  4.4  shows  extinction  of  the  solution  occurs 
at  t  =  1.305.  For  these  examples  we  have  chosen  Ax  =  0.05  and  At  =  0.001  because  this 
makes  our  stability  condition  <  1.25.  This  condition  is  always  true  for  cases  where 

^  =  vP  with  p  >  1,  along  with  our  initial  condition. 
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Figure  4.3.  Solution  to  equation  (1.23)  with  ^  =  u  and  F  =  ^/u.  Ax  =  0.05  and 
At  =  0.001. 

This  means  that,  as  before,  raising  the  power  of  one  variable  while  the  other  remains 
constant  will  increase  the  extinction  time  of  the  solution.  Although  the  solution  becomes 
extinct  at  a  later  time,  it  still  becomes  extinct  as  it  is  supposed  to.  We  see  here  that  ^ 
seems  to  have  a  much  greater  influence  on  the  solution  than  does  F. 

So,  if  increasing  the  power  of  or  F  in  (4.2)  increases  extinction  time,  then  making 
both  p  and  q  less  than  1  should  cause  extinction  to  occur  at  an  earlier  time.  From  Figure 
4.5  we  see  that  with  $  =  ^/u  and  F  =  ^/u  we  get  an  extinction  time  of  t  =  0.192. 

We  note  here  that  with  p  and  q  both  greater  than  1,  we  do  not  have  extinction  until 
infinity.  Although  we  do  not  show  any  plots  of  situations  such  as  this,  the  code  will  run 
until  the  solution  becomes  so  small  that  the  computer  rounds  the  solution  to  zero.  The 
program  that  was  used  to  run  this  code  starts  to  round  to  zero  at  10“^^^. 


4-7 
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Figure  4.4.  Solution  to  equation  (1.23)  with  and  F  =  y/u.  Ax  =  0.05  and 

At  =  0.001. 


This  is  not  promising.  The  fact  that  the  program  will  eventually  round  the  solution 
to  zero,  makes  it  hard  to  tell  whether  extinction  occurs  at  infinity  or  at  a  finite  time.  This 
especially  comes  into  play  when  we  have  a  and  F  that  have  a  large  finite  extinction  time. 

As  examples  of  solutions  with  infinite  extinction  times,  we  choose  =  tt  and  F  =  0, 
and  =  tt  and  F  =  u  because  their  solutions  reach  10“^^^  in  the  least  number  of  iterations 
and  therefore,  they  have  the  smallest  t  for  which  the  program  would  round  to  a  zero  solution 
of  all  the  solutions  that  become  extinct  at  infinity.  The  solutions  to  these  particular 
problems  are  rounded  to  zero  at  just  over  t  =  100. 

On  a  more  promising  note,  we  have  observed  that  when  finite  extinction  time  is 
supposed  to  occur,  the  observed  extinction  times  are  very  small  compared  to  t  =  100.  In 
a  previous  example,  our  numerical  solution  with  and  F  =  ^/u  goes  to  zero  at 

10“®.  This  means  that  for  the  iteration  just  before  the  numerical  solution  goes  to  zero. 
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Figure  4.5.  Solution  to  equation  (1.23)  with 
At  =  0.001. 


1 


the  numerical  solution  is  on  the  order  of  10  So,  the  numerical  solution  goes  to  zero  at 
a  much  larger  solution  than  10“^^^. 

Obviously  we  would  like  to  test  a  case  with  a  large  finite  extinction  time  such  as 
$  —  ^100  p  _  ^0.999  ggg  -^vhat  happens.  The  problem  with  this  is  that  the  code 
cannot  handle  such  a  case.  The  aforementioned  case  creates  ill-conditioned  matrices  for 
which  a  solution  cannot  be  found.  The  code  can  handle  and  F  =  From 

Figure  4.6  we  see  that  we  have  an  extinction  time  at  t  =  3.390.  Also,  the  iteration  just 
before  the  numerical  solution  goes  to  zero  is  on  the  order  of  10“^^.  Since  the  code  can 
only  handle  a  limited  number  of  cases,  and  the  limits  of  these  cases  have  extinction  times 
far  less  than  t  =  100,  then  we  feel  justified  in  saying  that  we  can  predict  whether  finite 
extinction  time  will  occur  for  any  given  case  that  the  code  can  handle. 
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3  ^><r 

3.5  0 

Figure  4.6.  Solution  to  equation  (1.23)  with  and  F  =  Ax  =  0.05  and 

At  =  0.001. 

The  code  can  handle  any  ^  =  u^,  where  p  >  ^.  It  can  also  handle  any  F  = 
where  q  >0.  The  only  exception  to  this  is  if  the  difference  between  p  and  q  is  large.  Also, 
remember  that  the  stability  condition  must  be  satisfied.  As  a  final  example,  Figure  4.7 
shows  and  F  =  u'^.  This  case  has  a  finite  extinction  time  at  t  =  0.413. 

4-3  Comparison  of  the  absorption-diffusion  problem  and  the  diffusion-free  problem 

In  this  section  we  will  compare  the  numerical  solution  of  the  one-dimensional  version 
of  equation  (1.23)  and  the  numerical  solution  of  (1.24)  with  M  =  1.  Equation  (1.23)  has 
the  same  boundary  and  initial  conditions  as  in  Section  4.2.  Since  we  are  concerned  with 
whether  the  solution,  u{x,t),  of  (1.23)  is  less  than  the  solution,  z(t),  of  (1.24),  then  we  will 
only  compare  the  maxj{u(xj,  t)}  to  z{t),  at  each  t.  In  this  thesis,  we  have  noticed  that  for 
the  equations  and  initial  conditions  considered,  the  maxj{tt(xi,  t)}  =  u{0.5,t). 
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Figure  4.7.  Solution  to  equation  (1.23)  with  and  F 

At  =  0.001. 


=  u 


Ax  =  0.1  and 
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Since  equation  (1.24)  has  finite  extinction  time  if  and  only  if  equation  (1.25)  holds, 
we  will  look  at  cases  where  and  F  satisfy  (1.25)  first.  For  the  sake  of  comparison,  we 
will  look  at  the  same  cases  as  we  did  in  Section  4.2. 

We  notice  that  in  Figures  4.8,  4.9,  4.10,  4.11,  4.12,  4.13,  and  4.14,  u  <  z.  This  is 
exactly  what  we  wanted,  and  is  a  good  sign  that  the  solution  of  (1.24)  becomes  extinct  in 
finite  time  when  (1.25)  holds. 

Now,  we  turn  to  cases  where  there  is  no  finite  extinction  time.  We  will  look  at  =  rt 
and  F  =  u,  and  and  F  =  u^.  Figure  4.15  shows  the  former.  Although  we  only 

show  the  first  3  time  units,  the  numerical  solution  does  not  get  rounded  to  zero  until  after 
t  =  100,  which  was  our  full  run-time.  We  found  that  u  <  z  for  all  t. 

Turning  to  the  case  where  and  F  =  we  clearly  see  from  Figure  4.16  that 

u  <  z.  This  solution  is  rounded  to  zero  well  after  our  run-time  of  t  =  100,  and  u  <  z  for 
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Figure  4.8.  Comparison  of  solutions  of  equation  (1.23)  and  (1.24)  with  ^  =  ^/u  and 
F  =  u.  Ax  =  .1  and  At  =  .001. 

the  entire  time.  This  means  that  we  have  found  u  <  z  for  every  case  that  we  tested,  no 
matter  when  the  extinction  time.  This  result  is  even  better  then  we  expected,  and  is  a 
good  sign  that  u  <  z  for  all  and  F  with  the  same  boundary  and  initial  conditions. 


V.  Conclusion 


5.1  Summary 

In  this  thesis  we  have  investigated  the  motivation  for  the  n-dimensional  absorption- 
diffusion  equation.  In  doing  so,  we  mathematically  developed  the  one-dimensional,  non¬ 
linear  absorption-diffusion  equation.  We  also  proved  analytically  that  for  every  absorption- 
diffusion  equation  with  $  and  F  satisfying  certain  conditions,  extinction  will  always  occur 
as  t  ^  oo. 

Next,  we  looked  at  various  examples  of  related  problems.  With  heat  pipe  problems, 
we  observed  the  development  of  a  computational  method  based  on  finite-difference  approx¬ 
imations.  We  then  considered  the  dead  core  problem,  which  is  a  real  world  example  of  finite 
extinction  time.  Finally,  we  reviewed  previous  work  that  was  done  on  finite  extinction  time 
problems  for  parabolic  equations. 

We  then  developed  a  numerical  method  in  order  to  calculate  solutions  to  the  one¬ 
dimensional,  non-linear  absorption-diffusion  equation.  We  tested  the  accuracy  of  this 
method  on  a  one-dimensional,  linear  diffusion  equation  with  a  solution  which  could  be 
written  in  closed  form,  and  we  saw  that  the  numerical  solution  generated  only  a  moderate 
amount  of  error.  With  such  a  small  amount  of  error  in  our  method,  we  were  confident 
that  our  numerical  scheme  was  accurate  and  proceeded  to  test  various  types  of  diffusion 
and  absorption  terms  in  order  to  see  which  ones  produced  extinction  in  finite  time.  Our 
extinction  results  were  consistent  with  the  analytical  results  of  Lair  and  Oxley  [10] . 

Next,  we  developed  a  numerical  method  to  computationally  solve  the  diffusion-free 
equation.  This  method  was  developed  in  a  similar  manner  to  the  method  for  the  one¬ 
dimensional,  non-linear  absorption-diffusion  equation.  This  ensured  that  we  had  a  high 
degree  of  accuracy  and  that  we  could  compare  the  numerical  solutions  of  the  two.  After 
comparison  of  the  numerical  solutions  with  our  particular  boundary  and  initial  conditions, 
we  found  that  for  all  of  the  cases  that  we  tested  and  no  matter  when  the  extinction  time, 
u  <  z. 

We  feel  confident  that  given  a  ‘h  and  F  for  which  we  can  carry  out  our  computational 
method,  we  can  predict  whether  or  not  finite  extinction  time  occurs.  Although  this  sounds 
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promising,  our  computational  method  is  limited  in  scope.  For  this  reason,  we  were  unable 
to  address  the  open  problem  where  we  have  slow  diffusion  and  weak  absorption,  yet  together 
they  are  strong. 

The  only  known  example  that  we  have  found  that  satisfies  (2.6)  is  where  ^  and  F  are 
piecewise  linear  functions  which  cross  each  other  infinitely  many  times  as  they  approach 
zero.  Therein  lies  the  problem,  in  that  a  computer  cannot  process  infinitely  many  pieces. 

5.2  Future  Work 

In  the  future,  we  would  like  to  increase  the  order  of  accuracy  of  our  approximations 
of  ut,  Ux,  and  Uxx  in  order  to  increase  the  accuracy  of  our  numerical  solutions.  We  would 
also  like  to  look  at  more  comparisons  between  u  and  z  in  order  to  further  justify  that  u  <  z 
in  all  cases  where  we  have  finite  extinction  time,  and  even  possibly  when  there  is  infinite 
extinction  time. 

We  would  also  like  to  modify  our  method  so  that  we  can  better  handle  cases  with 
fast  diffusion,  i.e.,  cases  where  and  do  not  exist  at  the  boundary.  By  doing  this,  we 
would  hope  that  we  could  make  our  Ax  increments  quite  a  bit  smaller  and  still  be  able  to 
obtain  a  numerical  solution,  which  would  be  more  accurate. 

Finally,  finding  a  different  example  of  and  F  which  satisfy  (2.6)  and  are  useable 
in  our  code  would  be  quite  helpful  in  providing  insight  into  the  open  problem.  Even  more 
useful  would  be  an  analytical  proof  of  whether  such  $  and  F  cause  finite  time  extinction. 
If  we  cannot  find  a  useable  ^  and  F,  it  would  be  beneficial  to  create  a  new  method  that 
can  handle  more  cases  with  finite  extinction  time,  in  order  to  do  a  more  in-depth  study  of 
fast  diffusion  and  strong  absorption  conditions. 
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Appendix  A.  Matlab  6.1  Code 


A.l  Absorption- Diffusion  Code 
function[z,i]  =heateqn(phi,fu,delx,delt,t) 
format  long 
syms  u 

phiprime=difF(phi,u); 
phidblprime=difF(phi,u,2) ; 

Nx=  [delx:delx:  1  — delx] ' ; 
indi=(l/delx)  — 1; 
indn=t/delt; 
z(:,l)=sin(pi*Nx); 

ppman=0; 

pdpman=0; 

fman=0; 

if  findsym(phiprime)=='  ' 

phiprime=double(phiprime) ; 
ppinit  =phiprime*ones(indi  ,indn+ 1 ) ; 
pp =phiprime  *ones  (indi  ,indn + 1 ) ; 
ppman=l; 

end 

if  findsym(phidblprime)=='  ' 

phidblprime=double(phidblprime); 

pdpinit=phidblprime*ones(indi,indn+l); 

pdp=phidblprime*ones(indi,indn+l); 

pdpman=l; 

end 

if  isa(fu,'double')==l 

funp=fu*ones  (indi  ,indn+ 1 ) ; 
fman=l; 

end 

hwin=waitbar(0, 'Please  wait...'); 
for  i=l:indn 

if  ppman==0 
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ppinit  ( :  ,i)  =subs(phiprime,  'u'  ,z( :  ,i) ) ; 


end 

if  pdpman==0 

pdpinit  (:  ,i)  =subs(phidblprime,  'u'  ,z( :  ,i) ) ; 

end 

ziplusl(l:indi— l,i)=z(2:indi,i); 
ziplusl(indi,i)=0; 
ziminusl(l,i)=0; 
ziminusl(2:indi,i)=z(l:indi— l,i); 

g=l-((2*delt)/(delx^2))*ppinit(:,i); 

h=(delt/(delx^2))*ppinit(:,i)+(delt/(4*delx^2))*pdpinit(:,i).*ziplusl(:,i)- 
(delt  /  (4*delx^  2) )  *pdpinit  ( :  ,i)  .*ziminusl  ( :  ,i) ; 
p= (delt  /  (debc^  2 ) )  *ppinit  ( :  ,i) + (delt  /  (4*delx^  2 ) )  *pdpinit  ( :  ,i)  •  *ziminus  1  ( :  ,i) 
(delt  /  (4*delx^  2) )  *pdpinit  ( :  ,i)  .*ziplusl  ( :  ,i) ; 

Q=diag(g); 

Q=Q+diag(h(l:indi-l),l); 

Q=Q+diag(p(2:indi),-l); 

if  fman==0 

funp(:,i)=subs(fu,'u',z(:,i)); 

end 

z(:,i+l)=Q*z(:,i)-delt*funp(:,i); 

[row,col]=find(z(:,i+l)<0); 
if  length(row)>=l 
z(row,i+l)=0; 

end 

[row,col]=find(z(:,i+l)==0); 
if  length(row)==indi 
i 

break 

end 

if  ppman==0 

pp(:,i+l)=subs(phiprime,'u',z(:,i+l)); 

end 
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if  pdpman==0 

pdp(:,i+l)=subs(phidblprime,'u',z(:,i+l)); 

end 

zp(l:indi— l,i+l)=z(2:indi,i+l); 

zp(indi,i+l)=0; 

zm(l,i+l)=0; 

zm(2:indi,i+l)=z(l:indi— l,i+l); 
a=l+(delt/(delx^2))*pp(:,i+l); 

b=(-delt/(2*delx^2))*pp(:,i+l)-(delt/(8*delx^2))*pdp(:,i+l).*zp(:,i+l)+  • 

(delt/(8*delx^2))*pdp(:,i+l).*zm(:,i+l); 
c=(-delt/(2*delx^2))*pp(:,i+l)-(delt/(8*delx^2))*pdp(:,i+l).*zm(:,i+l)+  . 

(delt/(8*delx^2))*pdp(:,i+l).*zp(:,i+l); 

Tone=diag(a); 

Tone=Tone+diag(b(  1  :indi— 1 ) ,  1 ) ; 

Tone=Tone+diag(c(2:indi),— 1); 

if  ppman==0 

pp(:,i)=subs(phiprime,'u',z(:,i)); 

end 

if  pdpman==0 

pdp(:  ,i)  =subs(phidblprime,  'u'  ,z( :  ,i) ) ; 

end 

zp(l:indi— l,i)=z(2:indi,i); 

zp(indi,i)=0; 

zm(l,i)=0; 

zm(2:indi,i)=z(l:indi— l,i); 
d= 1  -  (delt  /  (debc^  2 ) )  *pp( :  ,i) ; 

e=(delt/(2*delx^2))*pp(:,i)+(delt/(8*delx^2))*pdp(:,i).*zp(:,i)-  •  •  • 

(delt/(8*delx^2))*pdp(:,i).*zm(:,i); 
f=(delt/(2*delx^2))*pp(:,i)+(delt/(8*delx^2))*pdp(:,i).*zm(:,i)-  . . . 

(delt/(8*delx^2))*pdp(:,i).*zp(:,i); 

Ttwo=diag(d); 

Ttwo=Ttwo+diag(e(l:indi— 

Ttwo=Ttwo+diag(f(2:indi),— 1); 


if  fman==0 

funp(:,i+l)=subs(fu,'u',z(:,i+l)); 

end 

Rev=inv(Tone) ; 

Per=Rev*Ttwo; 

fud=funp(:,i+l)+funp(:,i); 

z(:,i+l)=Per*z(:,i)-Rev*[(delt/2)*fud]; 

[row,col]=find(z(:,i+l)<0); 
if  length(row)>=l 
z(row,i+l)=0; 

end 

[row,col]=find(z(:,i+l)==0); 
if  length(row)==indi 
i 

break 

end 

waitbar(i/indn,hwin) 

end 

close(hwin) 

A. 2  Diffusion- Free  Code 

function[z]  =zeqn(phi,fu,delt,t) 

format  long 

syms  u 

indn=t/delt; 

z(l)=l; 

h=waitbar(0, 'Please  wait...'); 
for  i=l:indn 

z(i+l)=z(i)— delt*(subs(fu,'u',z(i))+subs(phi,'u',z(i))); 

z(i+l)=z(i)— (delt/2)*(subs(fu,'u',z(i))+subs(phi,'u',z(i)))—  . . . 
(delt/2)*(subs(fu,'u',z(i+l))+subs(phi,'u',z(i+l))); 

waitbar(i/indn,h) 
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end 

close(h) 
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